module fmm_mod
   use kind_parameter
   implicit integer(I-N) ! 指定从A到F以及I，K开头的变量为整型数
   private  ! All entities are now module-private by default
   public lpointer,assigni,fmmmain! Explicitly export public entities

contains

   subroutine lpointer(lp, ln, maxia, ia, n, nexp, ntylr, ncellmx,&
      levmx, ligw, lrgw, nwksz, iwksz, nfield,&
      l_n, l_x, l_y, l_node, l_dnorm,&
      l_bc, l_a, l_b, l_xmax,&
      l_xmin, l_ymax, l_ymin, l_ielem, l_itree,&
      l_level, l_loct, l_numt, l_ifath, l_lowlev,&
      l_maxl, l_levmx, l_nexp, l_ntylr, l_tolerance,&
      l_ncellmx, l_nleafmx, l_mxl, l_u, l_ax,&
      l_sb, l_sx, l_ligw, l_lrgw, l_igwk,&
      l_rgwk, l_nwksz, l_iwksz, l_rwork, l_iwork,&
      l_xfield, l_nfield, l_f)

      implicit real*8(a-h,o-z)
      integer::ia(maxia)
      lp= 1

      l_n = l_address(1,maxia,ia,lp,4,1)
      l_x = l_address(2,maxia,ia,lp,8,n*2)
      l_y = l_address( 3,maxia,ia,lp,8,n*2)
      l_node = l_address( 4,maxia,ia,lp,4,n*2)
      l_dnorm = l_address( 5,maxia,ia,lp,8,n*2)
      l_bc = l_address( 6,maxia,ia,lp,8,n*2)
      l_a = l_address( 7,maxia,ia,lp,16,(nexp+1)*ncellmx)
      l_b = l_address( 8,maxia,ia,lp,16,(ntylr+1)*ncellmx)
      l_xmax = l_address( 9,maxia,ia,lp,8,1)
      l_xmin = l_address(10,maxia,ia,lp,8,1)
      l_ymax = l_address(11,maxia,ia,lp,8,1)
      l_ymin = l_address(12,maxia,ia,lp,8,1)
      l_ielem = l_address(13,maxia,ia,lp,4,n)
      l_itree = l_address(14,maxia,ia,lp,4,ncellmx)
      l_level = l_address(15,maxia,ia,lp,4,levmx+1)
      l_loct = l_address(16,maxia,ia,lp,4,ncellmx)
      l_numt = l_address(17,maxia,ia,lp,4,ncellmx)
      l_ifath = l_address(18,maxia,ia,lp,4,ncellmx)
      l_lowlev = l_address(19,maxia,ia,lp,4,1)
      l_maxl = l_address(20,maxia,ia,lp,4,1)
      l_levmx = l_address(21,maxia,ia,lp,4,1)
      l_nexp = l_address(22,maxia,ia,lp,4,1)
      l_ntylr = l_address(23,maxia,ia,lp,4,1)
      l_tolerance = l_address(24,maxia,ia,lp,8,1)
      l_ncellmx = l_address(25,maxia,ia,lp,4,1)
      l_nleafmx = l_address(26,maxia,ia,lp,4,1)
      l_mxl = l_address(27,maxia,ia,lp,4,1)
      l_u = l_address(28,maxia,ia,lp,8,n)
      l_ax = l_address(29,maxia,ia,lp,8,n)
      l_sb = l_address(30,maxia,ia,lp,8,n)
      l_sx = l_address(31,maxia,ia,lp,8,n)
      l_ligw = l_address(32,maxia,ia,lp,4,1)
      l_lrgw = l_address(33,maxia,ia,lp,4,1)
      l_igwk = l_address(34,maxia,ia,lp,4,ligw)
      l_rgwk = l_address(35,maxia,ia,lp,8,lrgw)

      l_nwksz = l_address(36,maxia,ia,lp,4,1)
      l_iwksz = l_address(37,maxia,ia,lp,4,1)
      l_rwork = l_address(38,maxia,ia,lp,8,nwksz)
      l_iwork = l_address(39,maxia,ia,lp,4,iwksz)
      l_xfield = l_address(40,maxia,ia,lp,8,nfield*2)
      l_nfield = l_address(41,maxia,ia,lp,4,1)
      l_f = l_address(42,maxia,ia,lp,8,nfield)
   end subroutine lpointer

!----------------------------------------------------------------------------

   integer function l_address(ln,maxia,ia,lp,ibyte,length)
      integer::ln,maxia,lp,ibyte,length,&
         iu,inc
      integer::ia(maxia)
      ! ln: 次序
      ! maxia: 允许的最大 variables 数目
      ! am = 大数组, 存储 SLATEC GMRES solver 的变量
      ! ia = 数组, 存储 am 中变量的位置
      ! lp; 数据起始位置
      ! ibyte; 单个数据的 byte 数
      ! length; 数据长度

      l_address = lp
      ia(ln) = lp

      iu = 16 ! 按照 byte*16 block 存储
      ! -1的作用是, 当恰好需要 iu 的整数倍时, 不会浪费1block
      inc = (ibyte*length-1)/iu+1
      lp = lp+inc
      if(ln .gt. maxia) then
         write(*,*)'!Specified # of variables maxia ', maxia ,' is too small'
         stop
      endif
   end function l_address

!----------------------------------------------------------------------------
   subroutine assigni(i,ii)
      integer(i4):: i,ii
      ii = i
      return
   end subroutine assigni

   subroutine assignd(d,dd)
      real(dp):: d,dd
      dd = d
      return
   end subroutine assignd

!----------------------------------------------------------------------------
   subroutine fmmmain(maxa, maxia, am, ia, n,x,y,node,dnorm,bc,&
   &a,b, xmax,xmin,ymax,ymin,ielem,itree,level,loct,numt,&
   &ifath,lowlev,maxl,levmx,nexp,ntylr,tolerance,ncellmx,&
   &nleafmx,mxl,u,ax,nfield,xfield,f,sb,sx,igwk,rgwk,&
   &ligw,lrgw,nwksz,iwksz,rwork,iwork)

      implicit real*8(a-h,o-z)
      complex(kind=8):: am(maxa), a,b
      dimension ia(maxia),ja(1),a(0:nexp,ncellmx),b(0:ntylr,ncellmx),&
         x(2,n),y(2,n),node(2,n),dnorm(2,n),bc(2,n),&
         ielem(n),itree(ncellmx),level(0:levmx),loct(ncellmx),&
         numt(ncellmx),ifath(ncellmx), u(n),ax(n),sb(n),sx(n),&
         igwk(ligw),rgwk(lrgw),rwork(nwksz),iwork(iwksz),&
         xfield(2,nfield),f(nfield)
      external matvec, msolve

! Input parameters and prepare the BEM model
      call prep_model(n,x,y,node,bc,dnorm,xfield,nfield,maxl,levmx,&
         nexp,ntylr,tolerance,xmin,xmax,ymin,ymax)

! Generate the quad-tree structure for the elements
      call tree(n,x,xmax,xmin,ymax,ymin,ielem,itree,level,loct,numt,&
         ifath,lowlev,maxl,levmx,ncellmx,nleafmx,nwksz,iwork)

! Compute the right-hand-side vector b with the FMM
      call fmmbvector(n,x,y,node,dnorm,bc,u,ax,a,b,xmax,xmin,ymax,ymin,&
         ielem,itree,level,loct,numt,ifath,&
         nexp,ntylr,ncellmx,lowlev,maxl,rwork,iwork)

! Solve the BEM system of equations Ax=b with the fast multipole BEM

! Prepare parameters for calling the iterative solver GMRES
! (SLATEC GMRES solver is used, which is available at www.netlib.org.
! See the documentation for the SLATEC GMRES solver for more information
! about the following related parameters)

      nelt = 1
      isym = 0
      itol = 0
      tol = tolerance
      iunit = 3
      igwk(1) = mxl
      igwk(2) = mxl
      igwk(3) = 0
      igwk(4) = 1
      igwk(5) = 10
      do i=1,n
         ax(i) = 0.d0
      enddo

      write(*,*) 'Call Equation Solver GMRES ...'
      call dgmres(n,u,ax, nelt,ia,ja,am,isym, matvec,msolve,itol,tol,&
         itmax,iter,er,ierr,iunit,sb,sx,rgwk,lrgw,igwk,ligw,&
         rwork,iwork)

      write(3,*) ' Error indicator from GMRES:', ierr
      write(*,*) ' Error indicator from GMRES:', ierr

! Output the boundary solution
      do i=1,n
         u(ielem(i)) = ax(i)
      enddo

      write(3,*) ' Fast Multipole BEM Solution:'
      do i=1,n
         write(3,*) i, u(i)
         write(7,*) i, u(i)
      enddo

! Evaluate the field inside the domain and output the results
      call domain_field(nfield,xfield,f,n,x,y,bc,node,dnorm,u)
      return
   end
!----------------------------------------------------------------------------

   subroutine prep_model(n,x,y,node,bc,dnorm,xfield,nfield,maxl,&
   &levmx,nexp,ntylr,tolerance,&
   &xmin,xmax,ymin,ymax)
      implicit real*8(a-h,o-z)
      dimension x(2,*),y(2,*),node(2,*),bc(2,*),dnorm(2,*),xfield(2,*)
      write(*,2) n, maxl, levmx, nexp, tolerance
      write(3,2) n, maxl, levmx, nexp, tolerance
2     format(' Total number of elements =', I12&
         //' Max. number of elements in a leaf =', I12&
         //' Max. number of tree levels =', I12&
         //' Number of terms used in expansions =', I12&
         //' Tolerance for convergence =', D12.3)
      write(*,*)
      write(3,*)

! Input the mesh data
      read(5,*)
      do i=1,n
         read(5,*) itemp, y(1,i), y(2,i)
      enddo
      read(5,*)
      do i=1,n
         read(5,*) itemp, node(1,i), node(2,i), bc(1,i), bc(2,i)
      enddo

! Input the field points inside the domain
      if (nfield .gt. 0) then
         read(5,*)
         do i=1,nfield
            read(5,*) itemp, xfield(1,i), xfield(2,i)
         enddo
      endif

! Compute mid-nodes and normals of the elements
      do i=1,n
         x(1,i) = (y(1,node(1,i)) + y(1,node(2,i)))*0.5
         x(2,i) = (y(2,node(1,i)) + y(2,node(2,i)))*0.5
         h1 = y(2,node(2,i)) - y(2,node(1,i))
         h2 = -y(1,node(2,i)) + y(1,node(1,i))
         el = sqrt(h1**2 + h2**2)
         dnorm(1,i) = h1/el
         dnorm(2,i) = h2/el
      enddo

! Determine the square bounding the problem domain (Largest cell used in FMM)
      xmin=x(1,1)
      xmax=x(1,1)
      ymin=x(2,1)
      ymax=x(2,1)

      do 10 i=2,n
         if(x(1,i).le.xmin) then
            xmin=x(1,i)
         elseif(x(1,i).ge.xmax) then
            xmax=x(1,i)
         endif
         if(x(2,i).le.ymin) then
            ymin=x(2,i)
         elseif(x(2,i).ge.ymax) then
            ymax=x(2,i)
         endif
10    continue

      scale = 1.05d0 ! Make the square slightly larger
      xyd = max(xmax-xmin,ymax-ymin)/2.d0
      xyd = xyd*scale
      cx = (xmin+xmax)/2.d0
      cy = (ymin+ymax)/2.d0
      xmin = cx-xyd
      xmax = cx+xyd
      ymin = cy-xyd
      ymax = cy+xyd

! Output nodal coordinates for plotting
      do i = 1,n
         write(8,*) x(1,i), x(2,i)
      enddo
      return
   end

!----------------------------------------------------------------------------
   subroutine bisec(x,ielem,n,xsep,nsep,ic)
      implicit real*8(a-h,o-z)
      dimension x(2,*),ielem(*)
      nsep = 1
      if(n.le.0) return

      do ifr=1,n
         if(x(ic,ielem(ifr)).le.xsep) then
            if(ifr.ne.nsep) then
               istore = ielem(nsep)
               ielem(nsep) = ielem(ifr)
               ielem(ifr) = istore
            endif
            nsep = nsep + 1
         endif
      enddo

      return
   end

!----------------------------------------------------------------------------
   subroutine tree(n,x,xmax,xmin,ymax,ymin,ielem,itree,level,loct,&
      numt,ifath,lowlev,maxl,levmx,ncellmx,nleafmx,&
      nwksz,iwork)
      implicit real*8(a-h,o-z)
      complex(dp):: a,b
      dimension x(2,*),ielem(*),itree(*),level(0:*),loct(*),numt(*),&
         ifath(*), iwork(*), nwk(4)

      do i=1,n
         ielem(i) = i ! Store the original element numbers in ielem
      enddo

! For the level 0 cell (largest cell)
      ifath(1) = 1 ! cell 的 parent cell
      itree(1) = 0 !
      loct(1) = 1 ! ith cell中, 单元的起始编号
      numt(1) = n ! number of elements
      level(0) = 1 ! level 起始编号
      level(1) = 2 !


      ndivx = 1
      lowlev = 1
      nleaf = 0
      nswa = 0

! For cells on level 1 to the lowest level (leaves)
      do 10 lev=1,levmx
         levp = lev-1
         levn = lev+1
         level(levn) = level(lev)

         if(level(lev).eq.level(levp)) goto 900

         ndivxp = ndivx
         ndivx = 2*ndivxp

         dxp = (xmax-xmin)/ndivxp ! Parent cell size
         dyp = (ymax-ymin)/ndivxp

         do 11 inp=level(levp),level(lev)-1
            itrp = itree(inp)

            if(numt(inp).gt.maxl.or.(lev.le.2 .and. numt(inp).ne.0))then
               itrpx = mod(itrp,ndivxp)
               itrpy = itrp/ndivxp
               xsep = xmin+(itrpx + 0.5d0)*dxp
               ysep = ymin+(itrpy + 0.5d0)*dyp

               call bisec(x,ielem(loct(inp)),numt(inp),ysep,nsepy, 2)
               call bisec(x,ielem(loct(inp)),nsepy-1,xsep,nsepx1,1)
               call bisec(x,ielem(loct(inp)+nsepy-1),numt(inp)-nsepy+1,xsep,nsepx2,1)

               nwk(1) = nsepx1-1
               nwk(2) = nsepy-nsepx1
               nwk(3) = nsepx2-1
               nwk(4) = numt(inp)-nsepy-nsepx2+2
               locc = loct(inp)

               do 12 icldy=0,1
                  do 12 icldx=0,1
                     icld = icldy*2+icldx+1

                     if(nwk(icld).gt.0) then
                        nrel = level(levn)
                        if(nrel.gt.ncellmx) then
                           write(*,*) " ncellmx error"
                           stop
                        endif
                        itree(nrel) = ((itrpy*2+icldy)*ndivxp + itrpx)*2 +icldx

                        loct(nrel) = locc
                        numt(nrel) = nwk(icld)
                        ifath(nrel) = inp
                        lowlev=lev

! Leaves:

                        if((lev.ne.1) .and.&
                        &(numt(nrel).le.maxl .or. lev.eq.levmx)) then
                           nleaf = nleaf+1
                           if(nleaf.gt.nleafmx) then
                              write(*,*) " nleafmx error"
                              stop
                           endif
                           nleaf3 = nleaf*3 - 1
                           iwork(nleaf3) = nrel ! Store cell number (icell)
                           iwork(nleaf3+1) = nswa + 1 ! Location of pre-cond'er
                           iwork(nleaf3+2) = 1 ! Initial value of switch isw
                           nswa = nswa + numt(nrel)**2
                           if(nswa.gt.nwksz) then
                              write(*,*) " nwksz error"
                              stop
                           endif
                        endif

                        level(levn) = nrel + 1
                        locc = locc + nwk(icld)
                     endif
12             continue
            endif
11       continue
10    continue

900   iwork(1) = nleaf ! Store number of leaves in iwork(1)
      write(3,*) ' Number of tree levels =', lowlev
      write(*,*) ' Number of tree levels =', lowlev
      write(3,*) ' Number of leaves =', nleaf
      write(*,*) ' Number of leaves =', nleaf
      write(3,*) ' Number of cells =', nrel
      write(*,*) ' Number of cells =', nrel
      write(3,*)
      write(*,*)
      return
   end



!----------------------------------------------------------------------------
   subroutine fmmbvector(n,x,y,node,dnorm,bc,u,ax,a,b,xmax,xmin,&
   &ymax,ymin,ielem,itree,level,loct,numt,ifath,&
   &nexp,ntylr,ncellmx,lowlev,maxl,rwork,iwork)
      implicit real*8(a-h,o-z)

      complex(dp):: a(0:nexp,ncellmx),b(0:ntylr,ncellmx)

      real(dp)::x(:,:),y(:,:),dnorm(:,:),bc(:,:),u(:),ax(:),rwork(:)

      integer:: node(:,:), ielem(:),itree(:),level(0:),loct(:),&
         numt(:),ifath(:),iwork(:)


! Switch the BC type
      do i=1,n
         if(bc(1,i) .eq. 1.) then
            bc(1,i) = 2.d0
         else
            bc(1,i) = 1.d0
         endif
      enddo

      do i=1,n
         u(i)= bc(2,ielem(i))
         ax(i) = 0.d0
      enddo

! Apply the FMM to conpute the right-hand side vector b
      call upward(u,n,y,node,&
         dnorm,bc,a,xmax,&
         xmin,ymax,ymin,ielem,&
         itree,level,loct,numt,&
         ifath,nexp,ncellmx,lowlev,maxl)

      call dwnwrd(u,ax,n,x,y,node,dnorm,bc,a,b,xmax,xmin,ymax,ymin,&
         ielem,itree,level,loct,numt,ifath,nexp,ntylr,ncellmx,&
         lowlev,maxl,rwork,iwork)

! Store b vector in u and switch the BC type back
      do i=1,n
         u(i) = - ax(i)
         if(bc(1,i) .eq. 1.) then
            bc(1,i) = 2.d0
         else
            bc(1,i) = 1.d0
         endif
      enddo
      return
   end

!----------------------------------------------------------------------------
   subroutine matvec(n,u,ax,nelt,ia,ja,am,isym)
      use iso_c_binding
      !implicit real*8(a-h,o-z)
      complex(dp):: a, b
      integer(i4)::n, ia(:), ja(:),isym
      real(dp):: u(:), ax(:)
      complex(dp):: am(:)

! Retrieve the pointers
      l_n = ia(1)
      l_x = ia(2)
      l_y = ia(3)
      l_node = ia(4)
      l_dnorm = ia(5)
      l_bc = ia(6)
      l_a = ia(7)
      l_b = ia(8)
      l_xmax = ia(9)
      l_xmin = ia(10)
      l_ymax = ia(11)
      l_ymin = ia(12)
      l_ielem = ia(13)
      l_itree = ia(14)
      l_level = ia(15)
      l_loct = ia(16)
      l_numt = ia(17)
      l_ifath = ia(18)
      l_lowlev = ia(19)
      l_maxl = ia(20)
      l_levmx = ia(21)
      l_nexp = ia(22)
      l_ntylr = ia(23)
      l_tolerance = ia(24)
      l_ncellmx = ia(25)
      l_nleafmx = ia(26)
      l_mxl = ia(27)
      l_u = ia(28)
      l_ax = ia(29)
      l_sb = ia(30)
      l_sx = ia(31)
      l_ligw = ia(32)
      l_lrgw = ia(33)
      l_igwk = ia(34)
      l_rgwk = ia(35)
      l_nwksz = ia(36)
      l_iwksz = ia(37)
      l_rwork = ia(38)
      l_iwork = ia(39)


! Evaluate matrix-vector multiplication Ax using the fast multipole BEM
!upward(u,n,y,node,dnorm
      call upward(am(l_n),am(l_y), am(l_node),&
         am(l_dnorm), am(l_bc), am(l_a), am(l_xmax),&
         am(l_xmin), am(l_ymax), am(l_ymin), am(l_ielem),&
         am(l_itree), am(l_level), am(l_loct), am(l_numt),&
         am(l_ifath), am(l_nexp), am(l_ncellmx), am(l_lowlev),&
         am(l_maxl))

      call dwnwrd(u,ax, am(l_n), am(l_x), am(l_y),&
         am(l_node), am(l_dnorm), am(l_bc), am(l_a),&
         am(l_b), am(l_xmax), am(l_xmin), am(l_ymax),&
         am(l_ymin), am(l_ielem), am(l_itree), am(l_level),&
         am(l_loct), am(l_numt), am(l_ifath), am(l_nexp),&
         am(l_ntylr), am(l_ncellmx), am(l_lowlev), am(l_maxl),&
         am(l_rwork), am(l_iwork))

      return
   end

!----------------------------------------------------------------------------
   subroutine msolve(n,r,z,nelt,ia,ja,am,isym,rwork,iwork)
      implicit real*8(a-h,o-z)
      complex*16 am(*)
      dimension r(*),z(*),ia(*),ja(1),rwork(*),iwork(*)

! Load the pointers
      l_loct = ia(16)
      l_numt = ia(17)

! Compute the preconditioning matrix
      call msolveinv(r,z,rwork,iwork, am(l_loct), am(l_numt))
      return
   end

!----------------------------------------------------------------------------
   subroutine msolveinv(r,z,rwork,iwork,loct,numt)
      implicit real*8(a-h,o-z)
      dimension r(*),z(*),iwork(*),rwork(*),loct(*),numt(*)

      nleaf = iwork(1)
      do l = 1,nleaf
         l3 = l*3-1
         inod = iwork(l3)

         indr = iwork(l3+1)
         indx = loct(inod)
         indi = indx+3*nleaf+1
         nr = numt(inod)

         ! https://netlib.org/lapack/explore-html/de/da4/group__double__blas__level1_ga21cdaae1732dea58194c279fca30126d.html
         ! DCOPY copies a vector, x, to a vector, y.
         ! dcopy (N,DX,INCX, DY,INCY)
         call dcopy(nr,r(indx),1,z(indx),1)
         call dluax(rwork(indr),nr,nr,z(indx),iwork(l3+2),&
         &iwork(indi),icon)

         if(icon.ne.0) then
            write(*,*) " dluax error, icon =", icon
            stop
         endif
         iwork(l3+2) = 2
      enddo

      return
   end

!----------------------------------------------------------------------------
! This subroutine solves linear system of equations Ax=b by LU decomposition.
!
! a .... given regular coefficient matrix.
! k .... given adjustable dimension for array a.
! n .... given order of matrix a.
! b .... given constant vector.
! isw .... given control information: c if 1, solve equations entirely.
! if 2, solve equations with last LU-decomposed entries.
! ip .... auxiliary 1 dimensioned array, size is n.
! transposition vector which represents
! row-exchanging by partial pivoting.
! icon.... resultant condition code. 
! if INFO = -i, the i-th argument had an illegal value

!
! Slave subroutines used (available at www.netlib.org):
! dgetrf, dgetrs
!
! dgetrf(M,N,A, LDA, IPIV, INFO)
! 
!    
!    
!    
!    
! 
! 
!    
!----------------------------------------------------------------------------

   subroutine dluax(a,k,n,b,isw,ip,icon)
      implicit real*8(a-h,o-z)
      dimension a(k,n),b(n),ip(n)
      data ione/1/

      icon = 30000
      if(isw.eq.1) go to 1000
      if(isw.eq.2) go to 1100
      go to 8000
1000  call dgetrf(n,n,a,k,ip,icon)
      isw = 2
1100  call dgetrs('n',n,ione,a,k,ip,b,n,icon)
8000  continue

      return
   end

!----------------------------------------------------------------------------
   subroutine upward(u,n,y,node,dnorm,&
      bc,a,xmax,xmin,ymax,ymin,ielem,&
      itree,level,loct,numt,ifath,&
      nexp,ncellmx,lowlev,maxl)

      complex(kind=8):: a(0:nexp,ncellmx),b, z0,zi
      integer(i4)::node(:,:)
      real(dp)::y(:,:),dnorm(:,:),bc(:,:),u(:),&
         xmax,xmin,ymax,ymin,&
         dx,dy,cx,cy,cxp,cyp
      integer(i4)::ielem(:),itree(:),level(0:),loct(:),numt(:),ifath(:)

      do i=1,level(lowlev+1)-1
         do k=0,nexp
            a(k,i) = (0.d0,0.d0) ! Clear multipole moments
         enddo
      enddo

      do 10 lev=lowlev,2,-1 ! Loop from leaf to level 2 cells (Upward)
         ndivx = 2**lev
         dx = (xmax-xmin)/ndivx ! Determine cell size
         dy = (ymax-ymin)/ndivx

         do 20 icell=level(lev),level(lev+1)-1 ! Loop for level l cells
            itr = itree(icell)
            itrx = mod(itr,ndivx)
            itry = itr/ndivx ! Position of the cell
            cx = xmin+(itrx + 0.5d0)*dx
            cy = ymin+(itry + 0.5d0)*dy ! Center of the cell

! Multipole expansion
            if(numt(icell) .le. maxl .or. lev .eq. lowlev) then ! Compute moment
               call moment(a(0,icell),y,node,ielem(loct(icell)),&
                  numt(icell),nexp,cx,cy,u(loct(icell)),&
                  bc,dnorm)
            endif

! M2M translation
            if(lev .ne. 2) then ! Do M2M translation to form moments
               cxp = xmin+(int(itrx/2)*2 + 1)*dx
               cyp = ymin+(int(itry/2)*2 + 1)*dy ! Center of parent cell
               z0 = cmplx(cx-cxp, cy-cyp,kind=dp) ! (z_c - z_c')
               io = ifath(icell) ! Cell no. of parent cell
               zi = (1.d0,0.d0)
               do k=0,nexp
                  do m=k,nexp
                     a(m,io) = a(m,io) + zi*a(m-k,icell) ! Use M2M
                  enddo
                  zi = zi*z0/(k+1)
               enddo
            endif

20       continue
10    continue
      return
   end

!----------------------------------------------------------------------------
   subroutine moment(a,y,node,ielem,num,nexp,cx,cy,u,bc,dnorm)
      complex(kind=8)::a(0:), z1,z2, zp1, zp2,zwbar,znorm
      integer(i4):: node(:,:),ielem(:)
      real(dp) :: q,phi,&
         y(:,:), u(:), bc(:,:),dnorm(:,:)

      do i=1,num ! Over elements in the leaf
         nelm = ielem(i) ! Element number
         n1 = node(1,nelm) ! Two ends of the element
         n2 = node(2,nelm)
         z1 = cmplx(y(1,n1)-cx, y(2,n1)-cy,kind=dp)
         z2 = cmplx(y(1,n2)-cx, y(2,n2)-cy,kind=dp)
         zwbar = conjg(z2 - z1)
         zwbar = zwbar/abs(zwbar) ! omega bar
         zp1 = z1*zwbar
         zp2 = z2*zwbar
         znorm = cmplx(dnorm(1,nelm),dnorm(2,nelm),kind=dp) ! complex normal n
         if(bc(1,nelm) .eq. 1.d0) then ! Assign values to phi and q
            phi = 0.D0
            q = u(i)
         else if(bc(1,nelm) .eq. 2.d0) then
            phi = u(i)
            q = 0.D0
         endif

! Compute moments:
         a(0) = a(0) - (zp2-zp1)*q ! G kernel
         do k=1,nexp
            a(k) = a(k) + (zp2-zp1)*znorm*phi ! F kernel
            zp1 = zp1*z1/(k+1)
            zp2 = zp2*z2/(k+1)
            a(k) = a(k) - (zp2-zp1)*q ! G kernel
         enddo
      enddo
      return
   end

!----------------------------------------------------------------------------
   subroutine dwnwrd(u,ax,n,x,y,node,dnorm,bc,a,b,xmax,xmin,&
   &ymax,ymin,ielem,itree,level,loct,numt,ifath,&
   &nexp,ntylr,ncellmx,lowlev,maxl,rwork,iwork)

      implicit real*8(a-h,o-z)
      complex*16 a,b, z0,zi,zo,zp

      dimension a(0:nexp,ncellmx),b(0:ntylr,ncellmx),x(2,*),y(2,*),&
      &node(2,*),dnorm(2,*),bc(2,*),u(*),ax(*),ielem(*),&
      &itree(*),level(0:*),loct(*),numt(*),ifath(*),&
      &rwork(*),iwork(*)

      data pi/3.141592653589793D0/
      pi2 = pi*2.d0
      do i=1,level(lowlev+1)-1
         do k=0,ntylr
            b(k,i) = (0.d0,0.d0)
         enddo
      enddo
      do i=1,n
         ax(i) = 0.d0
      enddo

      leaf = 0
      indr = 1
      indi = 1
      do 110 lev=2,lowlev ! Downward from level 2 cells to leaf cells
         ndivx = 2**lev
         dx = (xmax-xmin)/ndivx
         dy = (ymax-ymin)/ndivx
         do 120 icell=level(lev),level(lev+1)-1 ! Loop for level l cells
            itr = itree(icell)
            itrx = mod(itr,ndivx)
            itry = itr/ndivx ! Position of the cell
            cx = xmin+(itrx + 0.5d0)*dx
            cy = ymin+(itry + 0.5d0)*dy ! Center of the cell
            itrxp = itrx/2
            itryp = itry/2

! From the parent cell (use L2L)
            if(lev.ne.2) then
               cxp = xmin+(itrxp*2+1)*dx
               cyp = ymin+(itryp*2+1)*dy ! Center of the parent cell
               z0 = cmplx(cx-cxp, cy-cyp,kind=dp)
               io = ifath(icell) ! Cell no. of the parent cell
               zi = (1.d0,0.d0)
               do k=0,ntylr
                  do m=0,ntylr-k

                     b(m,icell) = b(m,icell) + zi*b(k+m,io) ! L2L translation
                  enddo
                  zi = zi*z0/(k+1)
               enddo
            endif

            do 130 jcell=level(lev),level(lev+1)-1
               jtr = itree(jcell)
               jtrx = mod(jtr,ndivx)
               jtry = jtr/ndivx
               jtrxp = jtrx/2
               jtryp = jtry/2

! The parents must be neighbours
               if(iabs(itrxp-jtrxp).gt.1 .or. iabs(itryp-jtryp).gt.1)&
               &goto 130

! For non-neighbours (cells in interaction list) (use M2L)
               if(iabs(itrx-jtrx).gt.1 .or. iabs(itry-jtry).gt.1) then
                  ccx = xmin + (jtrx + 0.5d0)*dx
                  ccy = ymin + (jtry + 0.5d0)*dy ! Center of the j cell
                  z0 = cmplx(cx-ccx, cy-ccy,kind=dp)
                  b(0,icell) = b(0,icell) - log(z0)*a(0,jcell)
                  zo = 1.
                  do m=1,nexp+ntylr
                     zo = zo/z0
                     kmin = max(0,m-nexp)
                     kmax = min(m,ntylr)
                     sgn = (-1.0)**kmin
                     do k=kmin,kmax
                        b(k,icell) = b(k,icell) + sgn*zo*a(m-k,jcell) ! M2L
                        sgn = -sgn
                     enddo
                     zo = zo*m
                  enddo

! Contribution from neighbouring leaves (use direct)
               elseif(numt(jcell).le.maxl .or.&
                  numt(icell).le.maxl .or. lev.eq.lowlev) then
                  if(icell.eq.jcell) then
                     leaf = leaf+1
                     leaf3 = leaf*3-1
                     if(iwork(leaf3).ne.icell) then
                        write(3,*) leaf,iwork(1),iwork(leaf3),icell,'check'
                        write(*,*) " icell error"
                        stop
                     endif
                     indr = iwork(leaf3+1)
                     indi = iwork(leaf3+2)
                  endif

                  call direct(ielem(loct(icell)),ielem(loct(jcell)),&
                  &node,x,y,numt(icell),numt(jcell),dnorm,&
                  &ax(loct(icell)),u(loct(jcell)),icell,jcell,&
                  &rwork(indr),indi,bc) ! Direct integration
               endif
130         continue

! Compute Ax if reach a leaf (Evaluate local expansion at each
! collocation point)
            if(numt(icell).le.maxl .or. lev.eq.lowlev) then
               fact = 1.d0
               do itylr=1,ntylr
                  fact = fact/itylr
                  b(itylr,icell) = b(itylr,icell)*fact
               enddo

               do in=1,numt(icell)
                  inax = loct(icell) + in-1 ! Element number in the tree
                  indx = ielem(inax) ! Original element number
                  zp = b(ntylr,icell)
                  z0 = cmplx(x(1,indx)-cx, x(2,indx)-cy,kind=dp)
                  do itylr=ntylr-1,0,-1
                     zp = zp*z0 + b(itylr,icell) ! Local expansion
                  enddo
                  zp = zp/pi2
                  ax(inax) = ax(inax) + dreal(zp) ! Array Ax
               enddo
            endif

120      continue
110   continue

      return
   end

!----------------------------------------------------------------------------
   subroutine direct(inod,jnod,node,x,y,ni,nj,dnorm,ax,u,icell,jcell,&
      amat,isw,bc)
      implicit real*8(a-h,o-z)
      dimension inod(*),jnod(*),node(2,*),x(2,*),y(2,*),&
      &dnorm(2,*),ax(*),u(*),amat(ni,*), bc(2,*)
      data pi/3.141592653589793D0/
      pi2 = pi*2.d0
      do j = 1, nj

         jind = jnod(j)
         al = sqrt((y(1,node(1,jind))- y(1,node(2,jind)))**2 +(y(2,node(1,jind)) &
            - y(2,node(2,jind)))**2) ! Element length

         do i = 1, ni
            iind = inod(i)
            x11 = y(1,node(1,jind)) - x(1,iind)
            x21 = y(2,node(1,jind)) - x(2,iind)
            x12 = y(1,node(2,jind)) - x(1,iind)
            x22 = y(2,node(2,jind)) - x(2,iind)
            r1 = sqrt(x11**2 + x21**2)
            r2 = sqrt(x12**2 + x22**2)
            d = x11*dnorm(1,jind) + x21*dnorm(2,jind)
            t1 = -x11*dnorm(2,jind) + x21*dnorm(1,jind)
            t2 = -x12*dnorm(2,jind) + x22*dnorm(1,jind)
            ds = abs(d)
            dtheta = datan2(ds*al,ds**2+t1*t2)
            aa = (-dtheta*ds + al + t1*log(r1) - t2*log(r2))/pi2
            if(d .lt. 0.d0) dtheta = -dtheta

            bb = -dtheta/pi2
            if(iind .eq. jind) bb = 0.5d0
            if(bc(1,jind) .eq. 1.) then
               ax(i) = ax(i) - aa*u(j)
               if(icell.eq.jcell .and. isw.eq.1) then ! Store coefficients in
                  amat(i,j) = - aa ! first iteration
               endif
            else if(bc(1,jind) .eq. 2.) then
               ax(i) = ax(i) + bb*u(j)
               if(icell.eq.jcell .and. isw.eq.1) then ! Store coefficients in
                  amat(i,j) = bb ! first iteration
               endif
            endif
         enddo
      enddo
      return
   end

!----------------------------------------------------------------------------
   subroutine domain_field(nfield,xfield,f,n,x,y,bc,node,dnorm,u)

      implicit real*8(a-h,o-z)
      dimension xfield(2,*), f(*), x(2,*),y(2,*),bc(2,*),node(2,*), dnorm(2,*),u(*)
      data pi/3.141592653589793D0/
      pi2 = pi*2.d0
      do i=1,nfield
         f(i) = 0.d0
      enddo
      do j=1,n ! Loop over all elements
         if(bc(1,j).eq.1) then
            f0 = bc(2,j)
            df0 = u(j)

         else if(bc(1,j).eq.2) then

            f0 = u(j)
            df0 = bc(2,j)
         endif
         al= sqrt((y(1,node(2,j)) - y(1,node(1,j)))**2 +&
            (y(2,node(2,j)) - y(2,node(1,j)))**2) ! Element length

         do i=1,nfield ! Loop over all field points inside the domain
            x11 = y(1,node(1,j)) - xfield(1,i)
            x21 = y(2,node(1,j)) - xfield(2,i)
            x12 = y(1,node(2,j)) - xfield(1,i)
            x22 = y(2,node(2,j)) - xfield(2,i)
            r1 = sqrt(x11**2 + x21**2)
            r2 = sqrt(x12**2 + x22**2)
            d = x11*dnorm(1,j) + x21*dnorm(2,j)
            t1 = -x11*dnorm(2,j) + x21*dnorm(1,j)
            t2 = -x12*dnorm(2,j) + x22*dnorm(1,j)
            ds = abs(d)
            theta1 = datan2(t1,ds)
            theta2 = datan2(t2,ds)
            dtheta = theta2 - theta1

            aa = (-dtheta*ds + al + t1*log(r1) - t2*log(r2))/pi2
            if(d .lt. 0.d0) dtheta = -dtheta
            bb = -dtheta/pi2
            f(i) = f(i) + aa*df0 - bb*f0
         enddo
      enddo

! Output results
      do i=1,nfield
         write(9,20) xfield(1,i), f(i)
      enddo
20    format(1x, 4E18.8)
      return
   end

!----------------------------------------------------------------------------
end module fmm_mod
